# import required modules
import pickle5
import pandas as pd
import numpy as np
DataPath = "../" # the folder containing the required files
OptFN = "{}options_full.ini".format(DataPath)
#### FUNCTIONS ####
# get run parameters from options file
def get_Options(OptFN, OptionsDic, AllowedPars):
try:
optionsfile=open(OptFN, "r")
for line in optionsfile :
if line[:1] != "#" and len(line.strip("\n".strip(""))) > 1 :
splitted_line = line.split(":")
key=splitted_line[0].split(" ")[0]
value=splitted_line[1].strip("\n").split(" ")[0]
if key in AllowedPars :
OptionsDic[key] = value
except:
print('Error reading %s file !' % (OptFN))
return OptionsDic
#retrieve list of lists based on 2 separators
def double_split(Text, sep1, sep2):
# further parse options
WholeL = []
MainL = Text.split(sep1)
for text in MainL:
WholeL.append(text.strip().split(sep2))
return WholeL
## load required data
# the options
OptionsDic = {}
AllowedPars = ["Partners"]
OptionsDic = get_Options(OptFN, OptionsDic, AllowedPars)
ComplexesL = double_split(OptionsDic["Partners"], ";", "~")[0]
# the AAvarsDic
AAvarsDic = {}
GlobalAAvarsDic = {}
with open(r"{}SR_AAvars_{}_vs_{}.pkl".format(DataPath,ComplexesL[0], ComplexesL[1]), "rb") as input_file:
GlobalAAvarsDic = pickle5.load(input_file)
## Drop irrelevant information
for complex in GlobalAAvarsDic.keys():
AAvarsDic = GlobalAAvarsDic[complex]
for key in AAvarsDic.keys():
del AAvarsDic[key]['IndBCData']
print("{} single mutant combinations found in SR_AAvars dictionary for complex partners {} (A) and {} (B) !".format(len(AAvarsDic.keys()), ComplexesL[0], ComplexesL[1]))
10956 single mutant combinations found in SR_AAvars dictionary for complex partners VN1551_VHH2 (A) and VN1554_TNFa-2 (B) !
dMutL = list(AAvarsDic.keys()) # list of 2D variant combination
# remove double mutants supported by less than Min_dBCs and SD more than NormEnrichFactor
Min_dBCs = 2
MaxSD_NormEnrichFactor = 0.8
for key in dMutL:
if ( AAvarsDic[key]['Stats'][2] >= Min_dBCs ) and ( AAvarsDic[key]['Stats'][1] <= (AAvarsDic[key]['Stats'][0] * MaxSD_NormEnrichFactor )):
AAvarsDic[key] = AAvarsDic[key]['Stats'] # Normalized enrichment, SD and diffent BCs supporting the combination of variants
else:
del AAvarsDic[key]
# print some information to see if everything is fine
i = 1
for key in AAvarsDic.keys():
if i <= 3:
print("{}:{}".format(key, AAvarsDic[key]))
i += 1
Y35A:Q88K:[0.3324051569820905, 0.07833378544510669, 4] N31E:L75G:[0.3874054823790009, 0.13775684010620373, 3] I57V:WT:[1.023368311213879, 0.14569355647951418, 149]
# evaluate number of mutants / partner and the maximum number of theoretical combinations
pA_Vars = set()
pB_Vars = set()
for key in AAvarsDic.keys():
key_split = key.split(":")
pA_Vars.add(key_split[0])
pB_Vars.add(key_split[1])
pA_VarsL = list(pA_Vars)
pA_PosL = [int(var[1:-1]) if var != "WT" else 0 for var in pA_VarsL]
pB_VarsL = list(pB_Vars)
pB_PosL = [int(var[1:-1]) if var != "WT" else 0 for var in pB_VarsL]
Tcomb = len(pA_VarsL) * len(pB_VarsL)
print("Total number of mutation in partner A: {}\nTotal number of mutation in partner B: {}\nMaximum theoretical number of 2D mutants: {}".format(len(pA_VarsL),len(pB_VarsL), Tcomb))
Total number of mutation in partner A: 122 Total number of mutation in partner B: 104 Maximum theoretical number of 2D mutants: 12688
# take position of an mutation
def takePos(muta):
if muta != "WT":
pos = int(muta[1:-1])
else:
pos = 0
return pos
pA_Vars = set()
pB_Vars = set()
for key in AAvarsDic.keys():
key_split = key.split(":")
pA_Vars.add(key_split[0])
pB_Vars.add(key_split[1])
pA_VarsL = PB_VarsL = []
pA_VarsL = list(pA_Vars)
pA_VarsL.sort(key=takePos)
pB_VarsL = list(pB_Vars)
pB_VarsL.sort(key=takePos)
# construct the pandas dataframe scaffold
AAvarsDF = pd.DataFrame(index=pA_VarsL, columns=pB_VarsL)
AAvarsSdDF = pd.DataFrame(index=pA_VarsL, columns=pB_VarsL)
AAvarsNdBCsDF = pd.DataFrame(index=pA_VarsL, columns=pB_VarsL)
AAvarsDF.shape
(122, 104)
# fill pandas dataframes (DF) with values
#enrichment DF
for key in AAvarsDic.keys():
key_split = key.split(":")
AAvarsDF.loc[key_split[0]][key_split[1]] = AAvarsDic[key][0]
#Standard deviation DF
for key in AAvarsDic.keys():
key_split = key.split(":")
AAvarsSdDF.loc[key_split[0]][key_split[1]] = AAvarsDic[key][1]
#Number of barcodes DF
for key in AAvarsDic.keys():
key_split = key.split(":")
AAvarsNdBCsDF.loc[key_split[0]][key_split[1]] = AAvarsDic[key][2]
AAvarsSdDF.head(10)
| WT | R32D | R32F | R32K | R32N | L75G | L75F | L75W | L75I | L75M | ... | E146L | S147N | S147M | S147F | S147T | S147E | S147D | S147K | S147G | S147L | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| WT | 0.280834 | 0.360511 | 0.23487 | 0.276151 | 0.327855 | 0.168185 | 0.415444 | 0.307319 | 0.267671 | 0.321664 | ... | 0.147814 | 0.185272 | 0.153633 | 0.235486 | 0.154807 | 0.387532 | 0.468021 | 0.206102 | 0.25404 | 0.33494 |
| S30R | 0.186665 | NaN | 0.03021 | 0.109499 | NaN | 0.181379 | 0.458014 | 0.144442 | 0.240737 | 0.295278 | ... | 0.109202 | 0.394844 | 0.396509 | 0.055291 | 0.617092 | NaN | 0.058589 | 0.276611 | 0.42216 | 0.036833 |
| S30L | 0.231614 | 0.191243 | NaN | 0.420196 | NaN | 0.020769 | 0.294963 | NaN | NaN | 0.108375 | ... | 0.117126 | 0.042798 | NaN | 0.234338 | NaN | NaN | 0.165077 | 0.025679 | 0.033494 | 0.021399 |
| S30K | 0.247644 | NaN | 0.027513 | 0.133331 | NaN | 0.330154 | NaN | NaN | NaN | 0.890193 | ... | 0.002093 | NaN | NaN | 0.608283 | NaN | NaN | NaN | 0.962949 | NaN | 0.285087 |
| S30G | 0.231018 | 0.601843 | 0.053333 | 0.148316 | NaN | 0.097118 | 0.303306 | 0.418674 | 0.19259 | 0.780512 | ... | 0.046539 | 0.053497 | 0.217199 | 0.091143 | 0.158752 | 0.528246 | NaN | 0.283485 | 0.20935 | 0.03824 |
| S30D | 0.160955 | NaN | 0.316432 | NaN | 0.519993 | 0.258905 | 0.622707 | NaN | NaN | 0.38952 | ... | 0.055077 | NaN | NaN | 0.0 | NaN | 0.169479 | NaN | 0.998124 | 0.444703 | 0.336821 |
| S30E | 0.29543 | NaN | 0.17208 | 0.012839 | NaN | 0.166962 | 0.131125 | 0.513573 | NaN | 0.0 | ... | 0.199997 | 0.492525 | 0.148066 | 0.357667 | 0.18189 | 0.207569 | NaN | 0.029629 | 0.166453 | NaN |
| S30F | 0.205494 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 1.626314 | ... | 0.057316 | NaN | NaN | 0.0 | 0.077036 | 0.082539 | NaN | NaN | NaN | 0.207404 |
| S30N | 0.327179 | NaN | 0.435768 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | 0.048432 | NaN | NaN | 0.503102 | NaN | NaN | NaN | NaN | 0.201761 | 0.0777 |
| N31D | 0.150951 | 0.491104 | 0.038602 | 0.03646 | NaN | 0.110898 | 0.104206 | NaN | 0.256786 | 0.065692 | ... | 0.011594 | 0.364185 | 0.051073 | 0.147361 | 0.431493 | 0.363867 | 0.165527 | 0.059917 | 0.006294 | 0.060528 |
10 rows × 104 columns
# prepare useful data for plotting
Enrich = AAvarsDF.transpose().to_numpy()
Sd = AAvarsSdDF.transpose().to_numpy()
NdBCs = AAvarsNdBCsDF.transpose().to_numpy()
tAAvarsDF = AAvarsDF.transpose()
# graphics plotting
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=go_colours.sequential.thermal
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Color sequential minus
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale =[[0.0, '#0d0887'],
[0.1111111111111111,'#46039f'],
[0.2222222222222222, '#7201a8'],
[0.3333333333333333, '#9c179e'],
[0.4444444444444444, '#bd3786'],
[0.5555555555555556, '#d8576b'],
[0.6666666666666666, '#ed7953'],
[0.7777777777777778, '#fb9f3a'],
[0.8888888888888888, '#fdca26'],
[1.0,'#f0f921']],
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Color Sequential from red to blue
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "rgb(165,0,38)"],
[0.1111111111111111, "rgb(215,48,39)"],
[0.2222222222222222, "rgb(244,109,67)"],
[0.3333333333333333, "rgb(253,174,97)"],
[0.4444444444444444, "rgb(254,224,144)"],
[0.5555555555555556, "rgb(224,243,248)"],
[0.6666666666666666, "rgb(171,217,233)"],
[0.7777777777777778, "rgb(116,173,209)"],
[0.8888888888888888, "rgb(69,117,180)"],
[1.0, "rgb(49,54,149)"]]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Color not Sequential from red to blue
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "rgb(165,0,38)"],
[0.0333333333333333, "rgb(215,48,39)"],
[0.0666666666666666, "rgb(244,109,67)"],
[0.1333333333333333, "rgb(253,174,97)"],
[0.2222222222222222, "rgb(254,224,144)"],
[0.3333333333333333, "rgb(224,243,248)"],
[0.4444444444444444, "rgb(171,217,233)"],
[0.5555555555555555, "rgb(116,173,209)"],
[0.8888888888888888, "rgb(69,117,180)"],
[1.0, "rgb(49,54,149)"]]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Color not Sequential from blue to red
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "rgb(49,54,149)"],
[0.0333333333333333, "rgb(69,117,180)"],
[0.0666666666666666, "rgb(116,173,209)"],
[0.1333333333333333, "rgb(171,217,233)"],
[0.2222222222222222, "rgb(224,243,248)"],
[0.3333333333333333, "rgb(254,224,144)"],
[0.4444444444444444, "rgb(253,174,97)"],
[0.5555555555555555, "rgb(244,109,67)"],
[0.8888888888888888, "rgb(215,48,39)"],
[1.0, "rgb(165,0,38)"]]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Color not Sequential from blue to red
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "rgb(49,54,149)"],
[0.0222222222222222, "rgb(69,117,180)"],
[0.0444444444444444, "rgb(116,173,209)"],
[0.0888888888888888, "rgb(171,217,233)"],
[0.1111111111111111, "rgb(224,243,248)"],
[0.2222222222222222, "rgb(254,224,144)"],
[0.3333333333333333, "rgb(253,174,97)"],
[0.4444444444444444, "rgb(244,109,67)"],
[0.5555555555555555, "rgb(215,48,39)"],
[1.0, "rgb(165,0,38)"]]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "rgb(166,206,227)"],
[0.0222222222222222, "rgb(69,117,180)"],
[0.0444444444444444, "rgb(31,120,180)"],
[0.0888888888888888, "rgb(178,223,138)"],
[0.1111111111111111, "rgb(51,160,44)"],
[0.2222222222222222, "rgb(251,154,153)"],
[0.3333333333333333, "rgb(253,174,97)"],
[0.4444444444444444, "rgb(244,109,67)"],
[0.5555555555555555, "rgb(215,48,39)"],
[1.0, "rgb(227,26,28)"]]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# graphics plotting
# Derived from Raph colors
import plotly.graph_objects as go
import plotly.colors as go_colours
heatmap_trace = go.Heatmap(
z = tAAvarsDF,
x = tAAvarsDF.columns,
y = tAAvarsDF.index,
customdata = np.stack((Enrich, Sd, NdBCs), axis=-1),
hovertemplate=('<b>{} '.format(ComplexesL[0]) + '(%{x})' + '<BR>{} '.format(ComplexesL[1]) +
'(%{y})</b><BR><BR>Enrichment: %{customdata[0]:,.2f}' + '<BR>sd: %{customdata[1]:,.3f}' + '<BR>Nbr_dBCs: %{customdata[2]}' +
'<extra></extra>'),
hoverongaps = False,
colorscale=[[0.0, "white"],
[0.0222222222222222, '#3b19ff'],
[0.0444444444444444, '#3b19ff'],
[0.0888888888888888, '#6fa8dc'],
[0.1333333333333333, '#93c47d'],
[0.2000000000000000, '#ffe599'],
[0.2777777777777777, '#f1c232'],
[0.3333333333333333, '#e69138'],
[0.4666666666666666, '#d63636'],
[0.8888888888888888, '#69031a'],
[1, '#69031a']]
)
fig = go.Figure(heatmap_trace)
fig['layout'].update(title='Normalized Enrichment Avarage',
yaxis={"title": ComplexesL[1]},
xaxis={"title": ComplexesL[0]},
width=1200,
height=800,
autosize=False,
)
fig.update_layout(xaxis_showgrid=False, yaxis_showgrid=False)
fig['layout']['yaxis']['autorange'] = "reversed"
fig.update_xaxes(side="top")
fig.show()
# print mutant related general information
print(AAvarsDic["W33R:N92D"])
[0.5317265029348482, 0.11034887139806233, 13]
# prepare to print mutant barcode specific information
WholeAAvarsDic = GlobalAAvarsDic[complex]
# print mutant combination related barcodes and their normalized enrichment
print(WholeAAvarsDic["W33R:N92D"])
[0.5317265029348482, 0.11034887139806233, 13]
print(WholeAAvarsDic["W33R:WT"])
[0.12334383285318813, 0.04871345301994174, 133]
print(WholeAAvarsDic["WT:N92D"])
[0.9348320424947547, 0.19905498618382542, 145]